Il movente

Lo studio nasce dalla curiosità di una studentessa di toccare con mano alcune delle misure macroeconomiche italiane più importanti, spesso menzionate da economisti e media.

I dati

Tutti i dati analizzati nel seguente studio sono stati scaricati dal portale FRED, “Federal Reserve Economic Data”, un vasto archivio di dati gestito dalla Federal Reserve Bank of St. Louis.
I dati scelti sono tutti rappresentati da serie storiche, tutte a frequenza trimestrale e non stagionalmente adattate.
La finestra temporale prescelta è quella che va dal primo trimestre del 1996 al quarto trimestre del 2022 e i dati sono tutti riferiti al contesto economico italiano.

Gli strumenti

Gli strumenti statistici utilizzati sono quelli appresi durante il corso di Laboratorio di Statistica, incentrato sull’analisi delle serie storiche e sull’interpretazione dei relativi risultati.

I packages R necessari per l’analisi che si intende condurre sono i seguenti.

library(forecast)
library(quantmod)
library(GGally)
library(tidyverse)
library(dygraphs)

Pulizia dei dati

Valori mancanti (missing values) o anomali (outliers) potrebbero avere effetti sui modelli che si andranno a costruire e sui relativi risultati. Risulta pertanto indispensabile verificare l’eventuale presenza di tali valori.

Missing values

anyNA(gdp)

## [1] FALSE

Outliers

tsoutliers(gdp)

## $index
## [1] 95 96 98
## 
## $replacements
## [1] 399055.6 397351.1 386198.6

La serie storica oggetto di studio non presenta missing values, pertanto non sarà necessario gestirli con metodi di imputazione; tuttavia, presenta outliers. I valori anomali presenti nella serie sono frutto degli effetti economici dell’emergenza Covid-19.

La funzione tsclean() utilizza la decomposizione robusta STL per la stima dei missing values e la sostituzione degli outliers. Possiamo utilizzarla per creare una serie storica in cui viene applicata la sostituzione degli outliers, in modo da poter confrontare i risultati dei metodi che verranno applicati sulle due serie.

gdp_clean <- tsclean(gdp)

autoplot(gdp_clean, series="Clean") +
  autolayer(gdp, series="Original") +
  ggtitle("Original series VS Clean series") +
  xlab("Year") +
  ylab("Currency (€)") +
  guides(colour=guide_legend(title="Series"))

Stagionalità

Si prosegue con la ricerca di un’eventuale stagionalità all’interno della serie storica del GDP.

ggsubseriesplot(gdp) +
  ylab("Currency (€)") + 
  ggtitle("Seasonal subseries plot: Real Gross Domestic Product")

ggseasonplot(gdp, year.labels = TRUE, year.labels.left = TRUE) +
  ylab("Currency (€)") + 
  xlab("Quarter")

  ggtitle("Seasonal plot: Real Gross Domestic")

## $title
## [1] "Seasonal plot: Real Gross Domestic"
## 
## attr(,"class")
## [1] "labels"

Come si poteva prevedere, il GDP può essere influenzato dalla stagionalità in quanto le attività economiche possono variare in modo regolare e prevedibile a seconda delle stagioni: il GDP mostra un andamento crescente il primo e il terzo trimestre, al contrario di ciò che accade il secondo e il quarto. Tuttavia, sembra che negli ultimi anni (orientativamente dal 2014) questo andamente sia stato meno netto e gli anni 2020 e 2021 hanno avuto andamenti anomali per le stesse considerazioni sopracitate legate all’emergenza Covid-19.

Decomposizione

É possibile decomporre la serie storica GDP isolando le diverse sue componenti: trend-ciclo, stagionalità e residuale. Si procede con la decomposizione STL (Seasonal and Trend decomposition using Loess), in quando è la più robusta ai valori anomali.

stl_fit <- stl(ts(as.vector(gdp), start=1996, frequency = 4), s.window = "periodic", robust=TRUE)
autoplot(stl_fit)

La decomposizione applicata ai dati originali appare fortemente influenzata dai valori anomali, nonostante la decomposizione STL sia la più robusta. Pertato, si applica la decomposizione alla serie gdp_clean.

stl_fit2 <- stl(ts(as.vector(gdp_clean), start=1996, frequency = 4), s.window = "periodic", robust=TRUE)
autoplot(stl_fit2)

La componente principale risulta essere quella di trend.

autoplot(gdp) +
  autolayer(stl_fit$time.series[,2], series="original data") +
  autolayer(stl_fit2$time.series[,2], series="clean data") +
  ylab("Currency (€)") +
  ggtitle("STL decomposition") +
  guides(colour=guide_legend(title="STL on"))

Possiamo, ora che sono state separate le componenti di trend e stagionalità, calcolarne la forza.

var_R <- var(stl_fit2$time.series[,3])
var_RT <- var(stl_fit2$time.series[,3] + stl_fit2$time.series[,2])
var_RS <- var(stl_fit2$time.series[,3] + stl_fit2$time.series[,1])

Forza del trend:

\[ F_t=max(0;1-\frac{Var(R_t)}{Var(R_t+T_t)}) \]

max(0,1-(var_R/var_RT))

## [1] 0.9661618

Forza della stagionalità:

\[ F_s=max(0;1-\frac{Var(R_t)}{Var(R_t+S_t)}) \]

max(0,1-(var_R/var_RS))

## [1] 0.8836213

Entrambi gli indici sono compresi tra 0 e 1. Essendo entrambi prossimi a 1, \(F_t=0.97\) e \(F_s=0.88\), entrambe le componenti sono forti, ma il trend risulta essere la componente principale.

Autocorrelazione

Per misurare la forza del legame lineare tra i valori della serie storica su due istanti di tempo diversi, t e t-h, dove h è detto ritardo o lag, si procede con lo studio dell’autocorrelzione della serie del GDP italiano.
L’autocorrelazione è data da

\[ \rho(h)= \frac{\gamma(h)}{\gamma(0)} \] dove

\[ \gamma(h) = \frac{1}{T} \sum_{t=1}^{T-h} (x_t - \bar{x})(x_{t+h} - \bar{x}) \]

e \(\gamma(0)\) si riduce a

\[ \gamma(0) = \frac{1}{T} \sum_{t=1}^{T} (x_t - \bar{x})^2 \] ovvero alla nozione di varianza.

L’insieme di tutti i valori di \(\rho(h)\) è detto funzione di autocorrelazione. Si tratta di una funzione simmetrica intorno allo zero.

ggAcf(gdp, plot=FALSE)

## 
## Autocorrelations of series 'gdp', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.551  0.697  0.393  0.658  0.240  0.403  0.137  0.428  0.082  0.283 
##     11     12     13     14     15     16     17     18     19     20 
##  0.061  0.365  0.030  0.213 -0.045  0.236 -0.075  0.099 -0.137  0.141

ggAcf(gdp) 

Dai valori delle autocorrelazioni, calcolati per diversi ritardi, risulta evidente che le autocorrelazioni di ritardo \(h=2\) sono quelle dal valore maggiore, coerentemente con le considerazioni fatte precedentemente circa la stagionalità della serie. Ciò si riflette anche nel grafico che mostra regolarità nel comportamento delle autocorrelazioni.
I valori delle autocorrelazioni permettono di rilevare la presenza di trend, in quanto inizialmente sono tutte positive e decacono a zero lentamente, e del relativo cambio di pendenza, in quanto c’è un cambio di direzione per i valori delle autocorrelazioni.
Nonostante i valori delle autocorrelazioni nella prima parte del grafico ricadano al di fuori dell’intervallo di confidenza (fissato al 95%), mostrando delle autocorrelazioni statisticamente significative, queste diventano statisticamente nulle nella parte finale, per valori elevati di h.

Previsioni

I dati di una serie storica, in questo caso del GDP italiano, possono essere utilizzati per fare previsione sul valore futuro che questa variabile assumerà. Chiaramente l’intervallo di confidenza diventa via via più ampio all’allontanarci dall’ultima osservazione a disposizione, rendendo le previsioni a lungo termine più incerte.

La presenza di valori anomali influisce notevolmente sulle capacità predittive dei modelli applicati. I modelli verranno addestrati sulle osservazioni dal 1996 al 2018 e testati sulle osservazioni dal 2019 al 2022 della serie storica trasformata gdp_clean per ridurre l’effetto dell’emergenza Covid sui dati.

train <- ts(gdp_clean[1:92], start=1996, frequency=4)
test <- window(gdp_clean, start=2019, end=c(2022,4))

mean_fit <- meanf(train, 16)
naive_fit <- naive(train, 16)
snaive_fit <- snaive(train, h=16)
drift_fit <- rwf(train, 16, drift=TRUE)

autoplot(gdp_clean) + 
  autolayer(mean_fit, series="Mean", PI=FALSE) + 
  autolayer(naive_fit, series="Naïve", PI=FALSE) +
  autolayer(snaive_fit, series="Seasonal Naïve", PI=FALSE) + 
  autolayer(drift_fit, series="Drift", PI=FALSE) +
  xlab("Year") + 
  ylab("Currency (€)") +
  ggtitle("Real Gross Domestic Product") +
  guides(colour=guide_legend(title="Forecast"))

rbind(
"Average" = accuracy(mean_fit, test)[2,c(2,3,5,6)],
"Naïve" = accuracy(naive_fit, test)[2,c(2,3,5,6)],
"Seasonal Naïve" = accuracy(snaive_fit, test)[2,c(2,3,5,6)],
"Drift" = accuracy(drift_fit, test)[2,c(2,3,5,6)])

##                    RMSE      MAE     MAPE     MASE
## Average        13891.99 11374.26 2.863320 1.672807
## Naïve          21352.64 16987.33 4.393137 2.498318
## Seasonal Naïve 13754.56 10821.43 2.785163 1.591503
## Drift          25741.42 22480.21 5.771592 3.306153

gdp_nc <- ts(gdp[1:96], start=1996, frequency=4)

train_nc <- ts(gdp_nc[1:84], start=1996, frequency=4)
test_nc <- window(gdp_nc, start=2017)

mean_fit_nc <- meanf(train_nc, 12)
naive_fit_nc <- naive(train_nc, 12)
snaive_fit_nc <- snaive(train_nc, h=12)
drift_fit_nc <- rwf(train_nc, 12, drift=TRUE)

autoplot(gdp_nc) + 
  autolayer(mean_fit_nc, series="Mean", PI=FALSE) + 
  autolayer(naive_fit_nc, series="Naïve", PI=FALSE) +
  autolayer(snaive_fit_nc, series="Seasonal Naïve", PI=FALSE) + 
  autolayer(drift_fit_nc, series="Drift", PI=FALSE) +
  xlab("Year") + 
  ylab("Currency (€)") +
  ggtitle("Real Gross Domestic Product") +
  guides(colour=guide_legend(title="Forecast"))

rbind(
"Average" = accuracy(mean_fit_nc, test_nc)[2,c(2,3,5,6)],
"Naïve" = accuracy(naive_fit_nc, test_nc)[2,c(2,3,5,6)],
"Seasonal Naïve" = accuracy(snaive_fit_nc, test_nc)[2,c(2,3,5,6)],
"Drift" = accuracy(drift_fit_nc, test_nc)[2,c(2,3,5,6)])

##                     RMSE       MAE     MAPE      MASE
## Average        12157.667 10544.807 2.585153 1.5139682
## Naïve           7254.374  6161.831 1.528049 0.8846834
## Seasonal Naïve 10155.249  9703.442 2.403862 1.3931695
## Drift           6999.073  5140.116 1.289655 0.7379909

spline_fit_nc <- splinef(train_nc, 12)

autoplot(spline_fit_nc)

autoplot(gdp_nc) + 
  autolayer(mean_fit_nc, series="Mean", PI=FALSE) + 
  autolayer(naive_fit_nc, series="Naïve", PI=FALSE) +
  autolayer(snaive_fit_nc, series="Seasonal Naïve", PI=FALSE) + 
  autolayer(drift_fit_nc, series="Drift", PI=FALSE) +
  autolayer(spline_fit_nc, series="Spline", PI=FALSE) +
  xlab("Year") + 
  ylab("Currency (€)") +
  ggtitle("Real Gross Domestic Product") +
  guides(colour=guide_legend(title="Forecast"))

rbind(
"Average" = accuracy(mean_fit_nc, test_nc)[2,c(2,3,5,6)],
"Naïve" = accuracy(naive_fit_nc, test_nc)[2,c(2,3,5,6)],
"Seasonal Naïve" = accuracy(snaive_fit_nc, test_nc)[2,c(2,3,5,6)],
"Drift" = accuracy(drift_fit_nc, test_nc)[2,c(2,3,5,6)],
"Spline" = accuracy(spline_fit_nc, test_nc)[2,c(2,3,5,6)])

##                     RMSE       MAE     MAPE      MASE
## Average        12157.667 10544.807 2.585153 1.5139682
## Naïve           7254.374  6161.831 1.528049 0.8846834
## Seasonal Naïve 10155.249  9703.442 2.403862 1.3931695
## Drift           6999.073  5140.116 1.289655 0.7379909
## Spline          7830.295  6396.832 1.595106 0.9184236

checkresiduals(drift_fit_nc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Random walk with drift
## Q* = 488.21, df = 8, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 8

checkresiduals(snaive_fit_nc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 108.08, df = 8, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 8

Real Exports of Goods and Services

Il Real Exports of Goods and Services si riferisce al valore totale di beni e servizi che l’Italia esporta verso altri Paesi, adeguato all’inflazione. Si tratta di un indicatore della competitività internazionale di un Paese.
Per l’Italia, le esportazioni sono una parte significativa dell’economia, poiché il Paese è noto per le sue esportazioni di prodotti di alta qualità.
Nota: Ci si riferirà a questa serie storica sempre con la sigla EXP.

Fonte: International Monetary Fund, Real Exports of Goods and Services for Italy [NXRNSAXDCITQ], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/NXRNSAXDCITQ, May 11, 2023.

exp <- ts(getSymbols("NXRNSAXDCITQ", src="FRED", auto.assign = FALSE), start=1996, frequency=4)

dyRangeSelector(dygraph(exp, 
                        main="Real Exports of Goods and Services",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Real Imports of Goods and Services for Italy

Il Real Imports of Goods and Services si riferisce al valore totale di beni e servizi che l’Italia importa da altri Paesi, adeguato all’inflazione. Si tratta di un indicatore economico della dipendenza dell’Italia dalle importazioni.
Come molte altre economie moderne, l’Italia importa una vasta gamma di beni e servizi.
Nota: Ci si riferirà a questa serie storica sempre con la sigla IMP.

Fonte: International Monetary Fund, Real Imports of Goods and Services for Italy [NMRNSAXDCITQ], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/NMRNSAXDCITQ, May 11, 2023.

imp <- ts(getSymbols("NMRNSAXDCITQ", src="FRED", auto.assign = FALSE), start=1996, frequency=4) 

dyRangeSelector(dygraph(imp, 
                        main="Real Imports of Goods and Services for Italy",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Early Estimate of Quarterly ULC Indicators: Total Labor Productivity

L’ULC, acronimo di Unit Labor Cost, è un indicatore macroeconomico che misura il costo del lavoro necessario per produrre un’unità di prodotto. L’ULC si calcola come il rapporto tra il costo del lavoro e la produttività del lavoro.
Questo indicatore può fornire informazioni importanti sulla capacità produttiva dell’economia italiana e sulla sua competitività internazionale.
Nota: Ci si riferirà a questa serie storica sempre con la sigla LAB.

Fonte: Organization for Economic Co-operation and Development, Early Estimate of Quarterly ULC Indicators: Total Labor Productivity for Italy [ULQELP01ITQ661N], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/ULQELP01ITQ661N, May 11, 2023.

lab <- ts(getSymbols("ULQELP01ITQ661N", src="FRED", auto.assign = FALSE) %>%
  window(start="1996-01-01"), start=1996, frequency=4) 

dyRangeSelector(dygraph(lab, 
                        main="Early Estimate of Quarterly ULC Indicators: Total Labor Productivity",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Real Final Consumption Expenditure for Italy

Con Real Final Consumption Expenditure ci si riferisce al valore totale della spesa finale effettuata dai residenti italiani per l’acquisto di beni e servizi, adeguato all’inflazione.
La spesa finale delle famiglie è un importante indicatore della domanda interna: una spesa finale elevata può indicare un aumento del reddito disponibile e una maggiore fiducia dei consumatori nell’economia. Tuttavia, una spesa finale eccessivamente elevata può anche indicare una maggiore inflazione o un aumento del debito delle famiglie.
Nota: Ci si riferirà a questa serie storica sempre con la sigla CONS.

Fonte: International Monetary Fund, Real Final Consumption Expenditure for Italy [NCRNSAXDCITQ], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/NCRNSAXDCITQ, May 11, 2023.

cons <- ts(getSymbols("NCRNSAXDCITQ", src="FRED", auto.assign = FALSE), start=1996, frequency=4) 

dyRangeSelector(dygraph(cons, 
                        main="Real Final Consumption Expenditure for Italy",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Real General Government Final Consumption Expenditure for Italy

Il Real General Government Final Consumption Expenditure si riferisce al valore totale della spesa finale effettuata dal governo per l’acquisto di beni e servizi (ad esempio le spese per la difesa, l’istruzione, la salute, la cultura, ecc.), adeguato all’inflazione.
La spesa pubblica può influire sull’economia nazionale attraverso il finanziamento di progetti infrastrutturali e di sviluppo, la creazione di posti di lavoro e l’aumento della domanda di beni e servizi. Tuttavia, una spesa pubblica eccessivamente elevata può anche portare a un aumento del debito pubblico e ad un aumento dei tassi di interesse.
Nota: Ci si riferirà a questa serie storica sempre con la sigla CONSGG.

Fonte: International Monetary Fund, Real General Government Final Consumption Expenditure for Italy [NCGGRNSAXDCITQ], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/NCGGRNSAXDCITQ, May 11, 2023.

consGG <- ts(getSymbols("NCGGRNSAXDCITQ", src="FRED", auto.assign = FALSE), start=1996, frequency=4)

dyRangeSelector(dygraph(consGG, 
                        main="Real General Government Final Consumption Expenditure for Italy",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Consumer Price Index: Total All Items for Italy

Il Consumer Price Index è un indicatore economico che misura la variazione dei prezzi di un paniere di beni e servizi consumati dalle famiglie. In pratica, l’indice dei prezzi al consumo fornisce una misura dell’inflazione in Italia.
L’aumento dei prezzi al consumo può avere un impatto significativo sulla capacità di acquisto delle famiglie e sull’andamento dell’economia.
Nota: Ci si riferirà a questa serie storica sempre con la sigla CPI.

Fonte: Organization for Economic Co-operation and Development, Consumer Price Index: Total All Items for Italy [CPALTT01ITQ657N], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/CPALTT01ITQ657N, May 10, 2023.

cpi <- ts(getSymbols("CPALTT01ITQ657N", src="FRED", auto.assign = FALSE) %>%
  window(start="1996-01-01", end="2022-10-01"), start=1996, frequency=4)

dyRangeSelector(dygraph(cpi, 
                        main="Consumer Price Index: Total All Items for Italy") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Real Gross Capital Formation for Italy

Con Real Gross Capital Formation ci si riferisce al valore totale degli investimenti effettuati dalle imprese, sia pubbliche che private, in beni di capitali al netto della deprezzamento e adeguato all’inflazione.
L’investimento in beni di capitale è una componente importante della spesa nazionale, poiché consente di aumentare la produttività e la capacità produttiva dell’economia italiana.
Nota: Ci si riferirà a questa serie storica sempre con la sigla CAP.

Fonte: International Monetary Fund, Real Gross Capital Formation for Italy [NIRNSAXDCITQ], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/NIRNSAXDCITQ, May 11, 2023.

cap <- ts(getSymbols("NIRNSAXDCITQ", src="FRED", auto.assign = FALSE), start=1996, frequency=4)

dyRangeSelector(dygraph(cap, 
                        main="Real Gross Capital Formation for Italy",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Costruzione del dataset

data <- cbind(gdp, exp, imp, lab, cons, consGG, cpi, cap)
colnames(data) <- c("gdp", "exp", "imp", "lab", "cons", "consGG", "cpi", "cap")

Per evitare i problemi legati alla multicollinearità, bisogna assicurarsi che le variabili non siano altamente correlate tra loro.

ggpairs(as.data.frame(data))

Risulta opportuno escludere le variabili imp e consGG, rispettivamente altamente correlate alle variabili exp e cons.

data <- cbind(gdp, exp, lab, cons, cpi, cap)
colnames(data) <- c("gdp", "exp", "lab", "cons", "cpi", "cap")

ggpairs(as.data.frame(data))

Dagli scatterplot risulta evidente l’esistenza di una relazione lineare tra la variabile GDP e le covariate prese in considerazione, suggerendo di poter costruire un modello di regressione lineare.

autoplot(data, facets = TRUE)

Modello di regressione lineare multiplo

Il modello più semplice che possiamo considerare è il modello di regressione lineare, che assume la presenza di un legame lineare tra la variabile gdp e le covariate sopra presentate. Inoltre, si aggiungono le variabili trend e season, legate alla struttura della serie storica dipendente.

\[ gdp = \beta_0 + \beta_1 exp + \beta_2 lab + \beta_3 cons + \beta_4 cpi + \beta_5 cap + \beta_6 trend + \beta_7 season + \epsilon \]

mod <- tslm(gdp ~ exp + lab + cons + cpi + cap + trend + season, data=data)

autoplot(gdp, series="Data") +
  autolayer(fitted(mod), series="Fitted") +
  ylab("GDP (€)")

Il modello sembra adattarsi molto bene ai dati, ma è necessario analizzare le sue statistiche di sintesi.

summary(mod)

## 
## Call:
## tslm(formula = gdp ~ exp + lab + cons + cpi + cap + trend + season, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4612.2 -1128.3   -94.5  1330.7  3438.6 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.443e+04  1.187e+04   2.058 0.042259 *  
## exp          4.583e-01  4.093e-02  11.196  < 2e-16 ***
## lab          4.127e+02  1.559e+02   2.647 0.009452 ** 
## cons         7.741e-01  2.479e-02  31.222  < 2e-16 ***
## cpi         -4.871e+02  3.709e+02  -1.313 0.192145    
## cap          5.126e-01  3.528e-02  14.532  < 2e-16 ***
## trend       -1.196e+02  3.050e+01  -3.922 0.000163 ***
## season2      5.416e+03  5.898e+02   9.182 7.21e-15 ***
## season3      4.230e+03  6.477e+02   6.531 2.91e-09 ***
## season4      5.847e+03  7.474e+02   7.822 6.07e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1891 on 98 degrees of freedom
## Multiple R-squared:  0.991,  Adjusted R-squared:  0.9902 
## F-statistic:  1202 on 9 and 98 DF,  p-value: < 2.2e-16

Un dato particolarmente rivelante è che tutte le variabili, ad eccezione del cpi (indice dei prezzi al consumo) sono statisticamente significative. Ciò è confermato anche dal fatto che il p-value della F-statistic ci porta a rifiutare l’ipotesi nulla \(H_0\) di nullità di tutti i predittori e ad accettare l’ipotesi alternativa \(H_1\) di almeno un predittore non nullo.

Il modello così costruito presenta un elevatissimo indice \(R^2\) corretto, una misura della bontà di adattamento. L’indice \(R^2\) indica la percentuale di varianza nella variabile dipendente spiegata dalle variabili indipendenti, ma non tiene conto del numero di variabili indipendenti utilizzate nel modello, il che può portare a una sovrastima dell’efficacia del modello quando vengono utilizzate molte variabili indipendenti. L’indice \(R^2\) corretto aumenta solo se l’aggiunta di una variabile indipendente migliora effettivamente la capacità del modello di spiegare la variazione nella variabile dipendente.
L’indice \(R^2\) corretto, così come quello standard, ha un valore compreso tra 0 e 1: quanto più il valore è prossimo a 1, tanto migliore è l’adattamento del modello ai dati. In questo caso l’indice corretto \(R^2= 0.99\) mostra un ottimo adattamento del modello ai dati.

Tuttavia, è indispensabile effettuare un’analisi dei residui. Dalle statistiche sopra riportate possiamo dedurre che la distribuzione dei residui è leggermente asimmetrica positiva e presenta uno standard error, misura della varianza dei residui, molto elevato. Si deduce che il modello di regressione non è in grado di spiegare tutta la variabilità presente nella variabile dipendente.

checkresiduals(mod)

## 
##  Breusch-Godfrey test for serial correlation of order up to 13
## 
## data:  Residuals from Linear regression model
## LM test = 58.143, df = 13, p-value = 1.126e-07

Modello di regressione lineare a tratti

Si potrebbe ipotizzare che questo problema sia legato ai cambi di pendenza nella serie storica del GDP, dunque si procede con la definizione di un modello lineare a tratti.

dyRangeSelector(dygraph(gdp, 
                        main="Real Gross Domestic Product",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Si procede con l’utilizzo di due nodi, uno in corrispondenza del quarto trimestre del 2007 e l’altro del primo trimestre del 2014.

t <- time(gdp)
tau.2007 <- c(2007,4)
tau.2014 <- c(2014,1)
tau1 <- ts(pmax(0, t - tau.2007), start = 1996)
tau2 <- ts(pmax(0, t - tau.2014), start = 1996)
mod_pw <- tslm(gdp ~ t + tau1 + tau2 + exp + lab + cons + cpi + cap + season)
summary(mod_pw)

## 
## Call:
## tslm(formula = gdp ~ t + tau1 + tau2 + exp + lab + cons + cpi + 
##     cap + season)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4244.5  -970.6   -86.9  1152.4  3952.4 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.475e+06  3.712e+05   3.973 0.000137 ***
## t           -7.121e+02  1.826e+02  -3.901 0.000178 ***
## tau1         1.836e+02  2.393e+02   0.768 0.444650    
## tau2        -3.777e+02  1.476e+02  -2.558 0.012090 *  
## exp          5.526e-01  4.418e-02  12.509  < 2e-16 ***
## lab         -1.806e+01  1.806e+02  -0.100 0.920590    
## cons         7.951e-01  2.621e-02  30.332  < 2e-16 ***
## cpi         -2.896e+02  3.477e+02  -0.833 0.406889    
## cap          5.238e-01  3.361e-02  15.582  < 2e-16 ***
## season2      3.963e+05  2.535e+05   1.563 0.121256    
## season3      4.090e+03  6.027e+02   6.786 9.48e-10 ***
## season4      3.977e+05  2.536e+05   1.568 0.120129    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1756 on 96 degrees of freedom
## Multiple R-squared:  0.9924, Adjusted R-squared:  0.9915 
## F-statistic:  1142 on 11 and 96 DF,  p-value: < 2.2e-16

checkresiduals(mod_pw)

## 
##  Breusch-Godfrey test for serial correlation of order up to 15
## 
## data:  Residuals from Linear regression model
## LM test = 43.828, df = 15, p-value = 0.000117

Il modello presenta ancora i residui correlati, ma il p-value del test di Breusch-Godfrey è aumentato, pur restando sotto il livello \(\alpha=0.05\). Si potrebbe, pertanto, pensare di aggiungere altri due noti, consapevoli del rischio di overfitting.

dyRangeSelector(dygraph(gdp_clean, 
                        main="Real Gross Domestic Product",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

t <- time(gdp)
tau.2007 <- c(2007,4)
tau.2014 <- c(2014,1)
tau.2018 <- c(2018,4)
tau.2020 <- c(2020,3)
tau1 <- ts(pmax(0, t - tau.2007), start = 1996)
tau2 <- ts(pmax(0, t - tau.2014), start = 1996)
tau3 <- ts(pmax(0, t - tau.2018), start = 1996)
tau4 <- ts(pmax(0, t - tau.2020), start = 1996)
mod_pw2 <- tslm(gdp ~ t + tau1 + tau2 + tau3 + tau2 + exp + lab + cons + cpi + cap + season)
summary(mod_pw2)

## 
## Call:
## tslm(formula = gdp ~ t + tau1 + tau2 + tau3 + tau2 + exp + lab + 
##     cons + cpi + cap + season)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4252.4  -966.2   -77.6  1126.5  4025.0 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.498e+06  3.727e+05   4.020 0.000116 ***
## t           -7.238e+02  1.833e+02  -3.949 0.000151 ***
## tau1         3.008e+02  2.754e+02   1.092 0.277495    
## tau2        -6.924e+02  3.936e+02  -1.759 0.081773 .  
## tau3         2.251e+02  2.609e+02   0.863 0.390465    
## exp          5.533e-01  4.425e-02  12.505  < 2e-16 ***
## lab         -9.249e+00  1.812e+02  -0.051 0.959391    
## cons         7.882e-01  2.746e-02  28.707  < 2e-16 ***
## cpi         -3.299e+02  3.513e+02  -0.939 0.350055    
## cap          5.359e-01  3.646e-02  14.696  < 2e-16 ***
## season2      3.421e+05  2.615e+05   1.308 0.193975    
## season3      4.200e+03  6.169e+02   6.809 8.84e-10 ***
## season4      3.435e+05  2.616e+05   1.313 0.192409    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1759 on 95 degrees of freedom
## Multiple R-squared:  0.9925, Adjusted R-squared:  0.9915 
## F-statistic:  1044 on 12 and 95 DF,  p-value: < 2.2e-16

checkresiduals(mod_pw2)

## 
##  Breusch-Godfrey test for serial correlation of order up to 16
## 
## data:  Residuals from Linear regression model
## LM test = 43.329, df = 16, p-value = 0.0002492

Il modello presenta ancora i residui correlati, ma il p-value del test di Breusch-Godfrey è aumentato, pur restando sotto il livello \(\alpha=0.05\).

Bisogna fare un’importante precisazione: la correlazione nei residui non rende le previsioni distorte, quindi queste non sono sbagliate, rende solo gli intervalli di previsione più ampi e quindi la stima più incerta.

rbind("lineare multiplo" = CV(mod),
      "lineare a tratti v.1" = CV(mod_pw),
      "lineare a tratti v.2" = CV(mod_pw2))

##                           CV      AIC     AICc      BIC     AdjR2
## lineare multiplo     4203785 1641.201 1643.951 1670.705 0.9901983
## lineare a tratti v.1 3757857 1627.032 1630.904 1661.899 0.9915435
## lineare a tratti v.2 3799479 1628.189 1632.705 1665.739 0.9915209

Siccome in termini di accuratezza previsiva tutti e tre i modelli presentati offrono ottime performance, il modello prescelto è il più semplice tra i due modelli lineari a tratti.

df <- as.data.frame(data)

df[,"Residuals"]  <- as.numeric(residuals(mod_pw))
p1 <- ggplot(df, aes(x=exp, y=Residuals)) +
  geom_point()
p2 <- ggplot(df, aes(x=lab, y=Residuals)) +
  geom_point()
p3 <- ggplot(df, aes(x=cons, y=Residuals)) +
  geom_point()
p4 <- ggplot(df, aes(x=cpi, y=Residuals)) +
  geom_point()
gridExtra::grid.arrange(p1, p2, p3, p4, nrow=2)

Previsioni

Su test set

É possibile utilizzare il modello di regressione lineare multiplo costruito mod per fare previsioni e valutare la bontà previsiva del modello. Per i motivi di cui sopra, il modello viene addestrato sulle osservazioni dal 1996 al 2016 e testato sulle osservazioni dal 2017 al 2019 della serie storica originale gdp.

nc_train <- ts(data[1:84,], start=1996, frequency = 4)
nc_test <- as.data.frame(ts(data[85:96,], start=2017, frequency = 4))

mod_nc <- tslm(gdp ~ exp + lab + cons + cpi + cap + trend + season, data=nc_train)

fcast_mod <- forecast(mod_nc, nc_test)

autoplot(gdp_nc) + 
  autolayer(fcast_mod, series="Forecast", PI=FALSE) +
  xlab("Year") + 
  ylab("Currency (€)") +
  ggtitle("Real Gross Domestic Product")

accuracy(fcast_mod, test_nc)[2,c(2,3,5,6)]

##         RMSE          MAE         MAPE         MASE 
## 2736.1123039 2360.1509437    0.5810829    0.3388581

Sul futuro

Possiamo, inoltre, utilizzare il modello così costruito per fare previsione sui valori futuri di GDP, previsioni che non possiamo testare.

fcast_mod_new <- forecast(mod$fitted.values, h=12)

autoplot(fcast_mod_new) +
  ggtitle("Forecasting") +
  xlab("Year") +
  ylab("Currency (€)")